PROGRAM translate
USE DFLIB
USE parm_mod
IMPLICIT NONE
INTEGER,PARAMETER :: maxxc = 2000,slen=300,maxpts=3000,maxbc=50,maxriv=10
INTEGER :: luin,luo1,luo2,luo4,istat,type_rm,lugis
INTEGER :: lineno,ndata,ic,icm,itemp1,itemp2,j,icinc,ikept,i,jj,nriv,nk,k1,k,itmp
INTEGER,DIMENSION(maxxc) :: xctype,n_gis
INTEGER,DIMENSION(maxxc) :: outXS       ! change cross section index into output index for each river
INTEGER,DIMENSION(maxxc) :: npt,num_rr,js,je,jbankL,jbankR
INTEGER,DIMENSION(maxriv) :: istart,iend,icOutput
INTEGER,DIMENSION(maxxc) :: nlev           ! number of levees in a XS
INTEGER,DIMENSION(maxxc) :: rating_num,rating_type
INTEGER,DIMENSION(maxxc) :: XC_lid_num
INTEGER,DIMENSION(maxxc) :: ndead,npdead           ! number of ineffective areas and permament ineffective areas in a XS
INTEGER,DIMENSION(maxxc) :: nblock           ! number of blocked areas in a XS
INTEGER,DIMENSION(maxriv) :: nki     ! number of internal boundary conditions
INTEGER,DIMENSION(maxriv,maxbc) :: KI           ! boundary condition type
INTEGER,DIMENSION(maxriv,maxbc) :: NXI          ! the cross section associated with the boundary condition
INTEGER,DIMENSION(maxxc) :: xsKept

CHARACTER(len=80) :: infile,outfile1,outfile2,outfile4,gisfile
CHARACTER(LEN=slen) :: line,des,answer,xcname,sdum
CHARACTER(len=slen),DIMENSION(maxriv) :: rivDes,reaDes
CHARACTER(len=slen),DIMENSION(maxxc) :: xcid
CHARACTER(len=slen),DIMENSION(maxxc,2) :: xcDes
CHARACTER(len=1),DIMENSION(maxbc,maxxc) :: typeDead  ! ineffective area tpye (T or F)
CHARACTER(len=300),DIMENSION(maxriv,maxbc) :: bcName,bcDes
CHARACTER(len=300),DIMENSION(20) :: line_char

LOGICAL :: floodPlain,skipxs,keepXS,read_gisfile,addpoints,addLowerPoint
LOGICAL :: exists
LOGICAL,DIMENSION(maxxc) :: sediment_present

REAL(prec) :: temp1,xloc,x1,y1,z1,rm1,rm2,min_dis
REAL(prec),DIMENSION(maxxc) :: xl,xr,bankl,bankr,fkec1,fkec2,skew,fminor
REAL(prec),DIMENSION(maxxc) :: sediment_elevation
REAL(prec),DIMENSION(maxbc,maxxc) :: sect_rr,loc_rr
REAL(prec),DIMENSION(0:maxxc) :: xt,x,xtL,xtR
REAL(prec),DIMENSION(maxpts,maxxc) :: y,z,x_gis,y_gis
REAL(prec),DIMENSION(maxpts,maxxc) :: q_rating,z_rating
REAL(prec),DIMENSION(maxpts,maxxc) :: lid_sta,lid_ele_lo,lid_ele_hi
REAL(prec),DIMENSION(maxpts,maxxc) :: yTmp,zTmp
REAL(prec),DIMENSION(maxbc,maxxc) :: leveel,leveer,hlev  ! levee locations and elevation
REAL(prec),DIMENSION(maxbc,maxxc) :: deadl,deadr,hdead  ! ineffective area locations and elevation
REAL(prec),DIMENSION(maxbc,maxxc) :: pdeadl,pdeadr,hpdead  ! ineffective area locations and elevation
REAL(prec),DIMENSION(maxbc,maxxc) :: blockl,blockr,hblock  ! blocked locations and elevation
REAL(prec),DIMENSION(maxriv,maxbc) :: XTI             ! dischance of the the boundary from cross section bcIc
REAL(prec),DIMENSION(maxriv,10,maxbc) :: bcVal             ! dischance of the the boundary from cross section bcIc
REAL(prec) :: ymax,deltZ
!REAL(prec),DIMENSION(maxxc) :: wse

icInc=1
lineno=0
ic=0
xt=0.; x=0.
nlev=0; leveel=0.; leveer=0.; hlev=0.
ndead=0; deadl=0.; deadr=0.; hdead=0.
npdead=0; pdeadl=0.; pdeadr=0.; hpdead=0.
nblock=0; blockl=0.; blockr=0.; hblock=0.
skew=0.;
fminor = 0.;
sediment_present = .FALSE.
sediment_elevation = 0.
floodPlain=.true.
nki=0.
icOutput=0
xcdes=''
n_gis = 0

WRITE(*,'(a)') 'SRH-1D TRANSLATE 4.0'
WRITE(*,'(a)') '===================='
WRITE(*,'(a)') 'This program transforms a HEC-RAS V4.1 geometry file to SRH-1D 4.0 format'
WRITE(*,'(a)') 'The translation of several structure types are not supported'
WRITE(*,*)

!******************************************************************************************
!     open all files
!******************************************************************************************
IF (NARGS() == 1) THEN
    DO I=1,10
		WRITE (*,'("Enter the HEC-RAS geometry file name: ",$)')
		READ (*,'(A)') infile
		infile = ADJUSTL(infile)
		INQUIRE (file=infile,EXIST=exists)
		IF (.NOT.exists) THEN
			WRITE (*,'(2a/)') '  <<ERROR>> FILE NOT FOUND: ', infile
			CYCLE
			IF(i==10)THEN
				WRITE (*,'(2a/)') ' EXITING....'
				STOP
			END IF
		ELSE
		    EXIT
		END IF
	ENDDO

	WRITE (*,'("Enter the name of the SRH-1D geometry file that will be created: ",$)')
	READ (*,'(A)') outfile2
	outfile2 = ADJUSTL(outfile2)

ELSEIF (NARGS() == 3) THEN
    CALL GETARG(1,infile)
    infile = ADJUSTL(infile)
	INQUIRE (file=infile,EXIST=exists)
	IF (.NOT.exists) THEN
		WRITE (*,'(2a/)') '  <<ERROR>> FILE NOT FOUND: ', infile, ' EXITING....'
	    STOP
	END IF
    CALL GETARG(2,outfile2)
    outfile2 = ADJUSTL(outfile2)
	INQUIRE (file=outfile2,EXIST=exists)
	IF (exists) THEN
!		WRITE (*,'(2a/)') '  <<ERROR>> FILE ALREADY EXISTS: ', outfile2, ' EXITING....'
!	    STOP
	END IF
END IF

outfile1 = 'out.dat'
outfile4 = 'thalweg.dat'

luin = 1
OPEN(UNIT=luin,file=infile,status='old',IOSTAT=istat,ACTION='READ')
IF(istat/=0)THEN
    WRITE(*,*) 'ERROR opening file',infile 
    STOP
END IF

!OPEN(UNIT=2,file="wse.dat",status='old',IOSTAT=istat,ACTION='READ')
!IF(istat/=0)THEN
!    WRITE(*,*) 'ERROR opening file',infile 
!    STOP
!END IF

luo1 = 11
OPEN(UNIT=luo1,file=outfile1,STATUS='unknown',IOSTAT=istat,ACTION='WRITE')
IF(istat/=0)THEN
    WRITE(*,*) 'ERROR opening file',outfile1
    STOP
END IF
luo2 = 12
OPEN(UNIT=luo2,file=outfile2,STATUS='unknown',IOSTAT=istat,ACTION='WRITE')
IF(istat/=0)THEN
    WRITE(*,*) 'ERROR opening file',outfile2
    STOP
END IF
luo4 = 14
OPEN(UNIT=luo4,file=outfile4,STATUS='unknown',IOSTAT=istat,ACTION='WRITE')
IF(istat/=0)THEN
    WRITE(*,*) 'ERROR opening file',outfile4
    STOP
END IF

!******************************************************************************************
! temperay wse
!******************************************************************************************
!DO i=1,81
!    read(2,*) wse(i)
!ENDDO


!******************************************************************************************
! YTT record:read title lines
!******************************************************************************************
YTT: DO 
  CALL commchk(line,lineno,ndata,luin,luo1)
  WRITE (luo1,'(i5,3x,A)') lineno,line
  IF (line(1:8)=='GEOM TIT') THEN
      WRITE (luo2,'(a,1x,a)') 'YTT',line(12:)
      continue
  ELSE
      EXIT YTT
  ENDIF  
ENDDO YTT

jj=0.
MAIN_LOOP: DO
!   ******************************************************************************************
!   TYPE: skip useless records and read type
!   ******************************************************************************************
	TYPE_sec: DO 
	    CALL commchk(line,lineno,ndata,luin,luo1)
	    WRITE (luo1,'(i5,3x,A)') lineno,line
		if(line(1:6)=='END OF') EXIT MAIN_LOOP
		IF(line(1:12)=='RIVER REACH=') THEN
	        jj = jj+1  ! increment xc #
			istart(jj)=ic+1
			nk = 0
			if(jj>maxriv) then
			    write(*,*) "exceeded maximum number of rivers"
			    write(*,*) "please contact Sedimentation and River Hydraulics Group of Reclamation to increase number of rivers"
				stop
			endif
		    CALL readLine(line(13:),line_char(1:2),2)
            rivDes(jj)=line_char(1)
            reaDes(jj)=line_char(2)
		    CYCLE TYPE_sec
	    ELSEIF (line(1:14)=='TYPE RM LENGTH') THEN
			READ(line(24:),'(i)')type_rm
			IF(type_rm==1)THEN
				ic = ic+1  ! increment xc #
				iend(jj)=ic
				if(ic>maxxc) then
					write(*,*) "exceeded maximum number of cross sections"
					write(*,*) "please contact Reclamation Sedimentation and River Hydraulics Group"
					stop
				endif
				IF(ndata>8)THEN
					READ(line(24:),*) xcType(ic),xcId(ic),xl(ic),x(ic),xr(ic)
				ELSE
					READ(line(24:),*) xcType(ic),xcId(ic)
				END IF
				xt(ic) = xt(ic-1)+x(ic-1)
			ELSE
				nki(jj) = nki(jj) + 1  ! increment internal #
				nk = nki(jj)
				READ(line(24:),*) type_rm,xti(jj,nk)
			END IF
		    EXIT TYPE_sec
		ELSE
			CYCLE TYPE_sec
		ENDIF  
	ENDDO TYPE_sec

!   ******************************************************************************************
!   DESCRIPTION: read descrition from Begin d to end d, or skip
!   ******************************************************************************************
	CALL commchk(line,lineno,ndata,luin,luo1)
	IF (line(1:6)=='BEGIN ') THEN
		WRITE (luo1,'(i5,3x,A)') lineno,line
		j=0      ! line of descriptions
		DESCRIPTION: DO 
			CALL commchk(line,lineno,ndata,luin,luo1)
			WRITE (luo1,'(i5,3x,A)') lineno,line
			IF (line(1:6)=='END DE') THEN
				EXIT DESCRIPTION
			ELSEIF(type_rm==1)THEN
			    j=j+1; if(j>2) j=2
				xcDes(ic,j)=line
				cycle DESCRIPTION
			ELSE
				bcDes(jj,nk)=line   ! only allow one line
				cycle DESCRIPTION
			ENDIF  
		ENDDO DESCRIPTION
	ELSE
        ! no description is read
	    BACKSPACE(luin); lineno=lineno-1
	ENDIF
!   ******************************************************************************************
!   skip Node name
!   ******************************************************************************************
	CALL commchk(line,lineno,ndata,luin,luo1)
    IF (line(1:3) == 'NOD') THEN
		WRITE (luo1,'(i5,3x,A)') lineno,line	
	ELSE
        ! no Nold Last Edited Time is read
	    BACKSPACE(luin); lineno=lineno-1
	ENDIF

!   ******************************************************************************************
!   DESCRIPTION: skip lines XS GIS Cut
!   ******************************************************************************************
	CALL commchk(line,lineno,ndata,luin,luo1)
	IF (line(1:6)=='XS GIS') THEN
		WRITE (luo1,'(i5,3x,A)') lineno,line
		READ(LINE(17:),*) n_gis(ic)
		DO k=1,n_gis(ic),2
		    READ(luin,'(4(g16.0))') (x_gis(j,ic),y_gis(j,ic),j=k,k+1)
            lineno=lineno+1
		    WRITE(luo1,'(i5,3x,4(f16.4))') lineno,(x_gis(j,ic),y_gis(j,ic),j=k,k+1)
		END DO
		!lineno = lineno + CEILING(REAL(n_gis(ic))/2.)
	ELSE
		! no XS GIS Cut
		BACKSPACE(luin); lineno=lineno-1
	ENDIF  
!   ******************************************************************************************
!   skip Node Last Edited Time
!   ******************************************************************************************
	CALL commchk(line,lineno,ndata,luin,luo1)
    IF (line(1:3) == 'NOD') THEN
		WRITE (luo1,'(i5,3x,A)') lineno,line	
	ELSE
        ! no Nold Last Edited Time is read
	    BACKSPACE(luin); lineno=lineno-1
	ENDIF
	
	SELECT CASE(type_rm)
	CASE(1_iprec)   !  it is station/elevation data
!		******************************************************************************************
!		read station and elevation data
!		******************************************************************************************
		CALL commchk(line,lineno,ndata,luin,luo1)
	    WRITE (luo1,'(i5,3x,A)') lineno,line
		IF(line(1:9) /= '#STA/ELEV') THEN
            write(*,*) "error: in reading file"
			CALL data_err('#STA/ELEV',line,lineno,luo1)
 		ENDIF
		READ(LINE(11:),*) npt(ic)
		READ(luin,'(10(f8.0))') (y(j,ic),z(j,ic),j=1,npt(ic))
        DO k=1,npt(ic),5
            lineno=lineno+1
		    WRITE(luo1,'(i5,3x,10(f8.2))') lineno,(y(j,ic),z(j,ic),j=k,k+4)
		END DO
		! lineno = lineno + CEILING(REAL(npt(ic))/5.)
!		******************************************************************************************
!		read manning's coefficients
!		******************************************************************************************
		if(ic==83) then
		    continue
		endif
		CALL commchk(line,lineno,ndata,luin,luo1)
	    WRITE (luo1,'(i5,3x,A)') lineno,line
		IF(line(1:6) == '#MANN=') THEN
			READ(LINE(7:),*,IOSTAT=istat) num_rr(ic)
            if(istat /= 0 ) CALL data_err('#MANN=',line,lineno,luo1)
            if(num_rr(ic)>0) then
			    READ(luin,'(9(f8.0))') ((loc_rr(j,ic),sect_rr(j,ic),temp1),j=1,num_rr(ic))
                DO k=1,num_rr(ic),3
                    lineno=lineno+1
		            WRITE(luo1,'(i5,3x,3(f8.2,f8.4,f8.2))') lineno,((loc_rr(j,ic),sect_rr(j,ic),temp1),j=k,k+2)
		        END DO
			endif
		ELSE
            write(*,*) "error: in reading file"
			CALL data_err('#Mann=',line,lineno,luo1)
 		ENDIF
!		******************************************************************************************
!		add levee location and elevation
!       if there is no levee, skip
!		******************************************************************************************
		CALL commchk(line,lineno,ndata,luin,luo1)
		IF(line(1:6) == 'LEVEE=') THEN
			WRITE (luo1,'(i5,3x,A)') lineno,line
			hlev(:,ic) = -99999.
			READ(LINE(7:),*,IOSTAT=istat) ((leveel(j,ic),leveer(j,ic),hlev(j,ic)),j=1,2)
            if(istat /= 0 ) THEN
			    READ(LINE(7:),*,IOSTAT=istat)                                             &
&                  leveel(1,ic),leveer(1,ic),hlev(1,ic),leveel(2,ic),leveer(2,ic)
				if(istat /= 0 ) CALL data_err('LEVEE=',line,lineno,luo1)
            ENDIF

			IF(abs(leveel(1,ic)-leveer(1,ic))<TINY) THEN  ! no left levee
			    nlev(ic)=1
				leveel(2,ic)=leveer(2,ic)
				leveer(2,ic)=y(npt(ic),ic)
				IF(hlev(2,ic)<TINY)THEN
					CALL bc_interp(z(:,ic),y(:,ic),leveel(2,ic),hlev(2,ic),npt(ic),2,2)
				END IF
				leveel(1,ic)=leveel(2,ic)
				leveer(1,ic)=leveer(2,ic)
				hlev(1,ic)=hlev(2,ic)
			ELSEIF(abs(leveel(2,ic)-leveer(2,ic))<TINY) THEN  ! no right levee
			    nlev(ic)=1
				leveel(1,ic)=y(1,ic)
				IF(hlev(1,ic)<TINY)THEN
					CALL bc_interp(z(:,ic),y(:,ic),leveer(1,ic),hlev(1,ic),npt(ic),2,2)
				END IF
            ELSE   ! both right and left levees
			    nlev(ic)=2
				leveel(1,ic)=y(1,ic)
				leveel(2,ic)=leveer(2,ic)
				leveer(2,ic)=y(npt(ic),ic)
				IF(hlev(1,ic)<TINY)THEN
					CALL bc_interp(z(:,ic),y(:,ic),leveer(1,ic),hlev(1,ic),npt(ic),2,2)
				END IF
				IF(hlev(2,ic)<TINY)THEN
					CALL bc_interp(z(:,ic),y(:,ic),leveel(2,ic),hlev(2,ic),npt(ic),2,2)
				END IF
			ENDIF
			DO j = 1,nlev(ic)
				leveel(j,ic) = MIN(y(npt(ic),ic),leveel(j,ic))
				leveer(j,ic) = MAX(y(1,ic),leveer(j,ic))
			END DO
		ELSE
			! no levee is read
			BACKSPACE(luin); lineno=lineno-1
 		ENDIF

!		******************************************************************************************
!		read ineffective location and elevation
!       if there is no ineffective area, skip
!		******************************************************************************************
		CALL commchk(line,lineno,ndata,luin,luo1)
		IF(line(1:6) == '#XS IN') THEN
			WRITE (luo1,'(i5,3x,A)') lineno,line
			READ(LINE(11:),*,IOSTAT=istat) ndead(ic)

			!CALL commchk(line,lineno,ndata,luin,luo1)
!			ndead(ic) = CEILING(REAL(ndata)/3.0)

            if(istat /= 0 ) CALL data_err('#XS IN',line,lineno,luo1)
            READ(luin,'(9(f8.0))') ((deadl(j,ic),deadr(j,ic),hdead(j,ic)),j=1,ndead(ic))
		    lineno = lineno + CEILING(REAL(ndead(ic)/3.))
			!READ(line,'(9(f8.0))')                                                       &
            !&				((deadl(j,ic),deadr(j,ic),hdead(j,ic)),j=1,ndead(ic))
			!lineno = lineno + 1
			CALL commchk(line,lineno,ndata,luin,luo1)
			WRITE (luo1,'(i5,3x,A)') lineno,line
			IF(line(1:6) == 'PERMAN') THEN
				CALL commchk(line,lineno,ndata,luin,luo1)
				READ(LINE,*,IOSTAT=istat) (typeDead(j,ic),j=1,ndead(ic))
				if(istat /= 0 .or. ndata/=ndead(ic) ) CALL data_err('typeDead',line,lineno,luo1)
				WRITE (luo1,'(i5,3x,A)') lineno,line
			ELSE
			    CALL data_err('PERMAN',line,lineno,luo1)
            ENDIF
  
 			DO j=1,ndead(ic)
			    if(deadl(j,ic)<TINY .and. deadr(j,ic)<TINY) typeDead(j,ic)='N'   ! it is not a dead area
				if(deadl(j,ic)<TINY) deadl(j,ic)=y(1,ic)
				if(deadr(j,ic)<TINY) deadr(j,ic)=y(npt(ic),ic)
				if(deadr(j,ic)<deadl(j,ic)) deadr(j,ic)=deadl(j,ic)+0.01
			ENDDO

			! seperate ineffective area with permament ineffective area
			itemp1=0;itemp2=0
 			DO j=1,ndead(ic)
                 IF(typeDead(j,ic)=='F') THEN     ! inefficient area
                     itemp1=itemp1+1
					 deadl(itemp1,ic)=deadl(j,ic)
					 deadr(itemp1,ic)=deadr(j,ic)
					 hdead(itemp1,ic)=hdead(j,ic)
                 ELSEIF(typeDead(j,ic)=='T') THEN     ! permament inefficient area
                     itemp2=itemp2+1
					 pdeadl(itemp2,ic)=deadl(j,ic)
					 pdeadr(itemp2,ic)=deadr(j,ic)
					 hpdead(itemp2,ic)=hdead(j,ic)
                 ELSEIF(typeDead(j,ic)=='N') THEN     ! this dead area is not valid

				 ELSE
				     write(*,*) 'wrong type of  inefficient area'
					 pause
					 stop
				 ENDIF
			ENDDO
			ndead(ic)=itemp1;npdead(ic)=itemp2
		ELSE
			! no dead is read
			BACKSPACE(luin); lineno=lineno-1
 		ENDIF

!		******************************************************************************************
!		read blocked location and elevation
!       if there is no blocked, skip
!		******************************************************************************************
		CALL commchk(line,lineno,ndata,luin,luo1)
		IF(line(1:6) == '#BLOCK') THEN
			WRITE (luo1,'(i5,3x,A)') lineno,line
			READ(LINE(17:),*,IOSTAT=istat) nblock(ic)
            if(istat /= 0 ) CALL data_err('#BLOCK',line,lineno,luo1)
			READ(luin,'(9(f8.0))') ((blockl(j,ic),blockr(j,ic),hblock(j,ic)),j=1,nblock(ic))
			lineno = lineno + 1
            IF(nblock(ic)==2) THEN
				IF(blockl(1,ic)<TINY) blockl(1,ic)= y(1,ic)
				IF(blockr(2,ic)<TINY) blockr(2,ic)= y(npt(ic),ic)
				IF(blockr(1,ic)<TINY) THEN        ! no left blocked area					
					blockl(1,ic)= blockl(2,ic)
					blockr(1,ic)= blockr(2,ic)
					hblock(1,ic)= hblock(2,ic)
					nblock(ic)=1
				ELSEIF(blockl(2,ic)<TINY) THEN        ! no right blocked area					
					nblock(ic)=1
				ENDIF
			ENDIF
		ELSE
			! no blocked area is read
			BACKSPACE(luin); lineno=lineno-1
 		ENDIF	

		DO j=1,nblock(ic)
			if(hblock(j,ic)<TINY) hblock(j,ic)=maxval(z(1:npt(ic),ic))
		ENDDO	
!		******************************************************************************************
!       #XS Lid 		
!		******************************************************************************************
		CALL commchk(line,lineno,ndata,luin,luo1)
		IF(line(1:7) == '#XS LID') THEN  ! Cross section lid
	        WRITE (luo1,'(i5,3x,A)') lineno,line
		    READ(line(9:),*) XC_lid_num(ic)
		    IF(XC_lid_num(ic)>0)THEN
       		    READ(luin,'(9(f8.0))') (lid_sta(j,ic),lid_ele_lo(j,ic),lid_ele_hi(j,ic),j=1,XC_lid_num(ic))
       		    lineno = lineno + CEILING(REAL(XC_lid_num(ic))/3.)
       		END IF
		ELSE
			BACKSPACE(luin); lineno=lineno-1
 		ENDIF
!		******************************************************************************************
!		read bank station number
!		******************************************************************************************
		CALL commchk(line,lineno,ndata,luin,luo1)
	    WRITE (luo1,'(i5,3x,A)') lineno,line
		IF(line(1:6) /= 'BANK S') THEN
            write(*,*) "error: in reading file"
			CALL data_err('Bank Sta=',line,lineno,luo1)
 		ENDIF
		READ(LINE(10:),*) bankl(ic),bankr(ic)

!		******************************************************************************************
!		XS Rating Curve= 
!		******************************************************************************************
		CALL commchk(line,lineno,ndata,luin,luo1)
		IF(line(1:9) == 'XS RATING') THEN  ! forced rating curve
	        WRITE (luo1,'(i5,3x,A)') lineno,line
		    READ(line(17:),*) rating_num(ic)! ,rating_type(ic)
		    IF(rating_num(ic)>0)THEN
       		    READ(luin,'(10(f8.0))') (q_rating(j,ic),z_rating(j,ic),j=1,rating_num(ic))
       		    lineno = lineno + CEILING(REAL(rating_num(ic))/2.)
       		END IF
		ELSE
			BACKSPACE(luin); lineno=lineno-1
 		ENDIF
!		******************************************************************************************
!		skip Sediment Elevation= 
!		******************************************************************************************
		CALL commchk(line,lineno,ndata,luin,luo1)
		IF(line(1:12) == 'SEDIMENT ELE') THEN  ! update geometry so that minimum sediment elevation is met
	        WRITE (luo1,'(i5,3x,A)') lineno,line
            READ(line(20:),*) sediment_elevation(ic) 
            sediment_present(ic) = .TRUE.
            DO k=1,npt(ic)
                z(k,ic) = MAX(z(k,ic),sediment_elevation(ic)) ! take maximum of sediment elevation and original bed elevation
            END DO
		ELSE
			BACKSPACE(luin); lineno=lineno-1
        ENDIF
!		******************************************************************************************
!		Skip XS HTab records (Starting El and Incr= , Horizontal distribution, etc...)
!		******************************************************************************************
        htab: DO
		    CALL commchk(line,lineno,ndata,luin,luo1)
		    IF(line(1:7) == 'XS HTAB') THEN  ! forced rating curve
    	        WRITE (luo1,'(i5,3x,A)') lineno,line
            ELSE
			    BACKSPACE(luin); lineno=lineno-1
                EXIT htab
            ENDIF
        END DO htab
!		******************************************************************************************
!		read skew if present
!		******************************************************************************************
		CALL commchk(line,lineno,ndata,luin,luo1)
		IF(line(1:4) == 'SKEW') THEN
	        WRITE (luo1,'(i5,3x,A)') lineno,line
			READ(LINE(12:),*) skew(ic)
		ELSE
			BACKSPACE(luin); lineno=lineno-1
        ENDIF
!		******************************************************************************************
!		read minor losses
!		******************************************************************************************
		CALL commchk(line,lineno,ndata,luin,luo1)
		IF(line(1:5) == 'MINOR') THEN
    	    WRITE (luo1,'(i5,3x,A)') lineno,line
    		READ(LINE(24:),*) fminor(ic)
        ELSE
			BACKSPACE(luin); lineno=lineno-1            
 		ENDIF
!		******************************************************************************************
!		read local contraction and expension loss
!		******************************************************************************************
		CALL commchk(line,lineno,ndata,luin,luo1)
	    WRITE (luo1,'(i5,3x,A)') lineno,line
		IF(line(1:6) /= 'EXP/CN') THEN
            write(*,*) "error: in reading file"
			CALL data_err('Exp/Cntr=',line,lineno,luo1)
 		ENDIF
		READ(LINE(10:),*) fkec1(ic),fkec2(ic)
    CASE(3_iprec)     !  bridge
!        write(*,*) "unknown type of data"
!        write(*,*) "please contact Victor Huang to add internal boundary conditions for bridge."
        ic=ic
		cycle main_loop
    CASE(5_iprec)     !  radial gate
!		******************************************************************************************
!		read IW Pilot Flow=0
!		******************************************************************************************
		nxi(jj,nk)=ic
		ki(jj,nk)=8         ! radial gate
		CALL commchk(line,lineno,ndata,luin,luo1)
	    WRITE (luo1,'(i5,3x,A)') lineno,line
		IF(line(1:6) /= 'IW PIL') THEN
            write(*,*) "error: in reading file"
			CALL data_err('IW Pilot Flow=',line,lineno,luo1)
 		ENDIF
!		******************************************************************************************
!		read #Inline Weir SE=
!		******************************************************************************************
		CALL commchk(line,lineno,ndata,luin,luo1)
	    WRITE (luo1,'(i5,3x,A)') lineno,line
		IF(line(1:12) /= '#INLINE WEIR') THEN
            write(*,*) "error: in reading file"
			CALL data_err('#Inline Weir SE=',line,lineno,luo1)
 		ENDIF
		READ(LINE(17:),*) itemp1
		READ(luin,'(10(f8.0))') (temp1,j=1,2*itemp1)
		lineno = lineno + 1
!		******************************************************************************************
!		read IW Dist,WD,Coef,Skew,MaxSub,Min_El,Is_Ogee,SpillHt,DesHd
!		******************************************************************************************
		CALL commchk(line,lineno,ndata,luin,luo1)
	    WRITE (luo1,'(i5,3x,A)') lineno,line
		IF(line(1:6) /= 'IW DIS') THEN
            write(*,*) "error: in reading file"
			CALL data_err('IW Dist',line,lineno,luo1) 
 		ENDIF
		READ(luin,*) xti(jj,nk)
		lineno = lineno + 1
!		*********************************************************************************************
!		read IW Gate Name Wd,H,Inv,GCoef,Exp_T,Exp_O,Exp_H,Type,WCoef,Is_Ogee,SpillHt,DesHd,#Openings
!		*********************************************************************************************
		CALL commchk(line,lineno,ndata,luin,luo1)
	    WRITE (luo1,'(i5,3x,A)') lineno,line
		IF(line(1:6) /= 'IW GAT') THEN
			BACKSPACE(luin); lineno=lineno-1			
		ELSE
		    CALL commchk(line,lineno,ndata,luin,luo1)
	        WRITE (luo1,'(i5,3x,A)') lineno,line
		    CALL readLine(line,line_char(1:16),16)
		    bcName(jj,nk)=line_char(1)
    		
		    i = 2; j = 2
		    IF(line_char(i)/="")THEN
		        READ(line_char(i),*) bcVal(jj,j,nk)
		    ELSE
		        bcVal(jj,j,nk) = 0.
		    END IF
    		
		    i = 3; j = 3
		    IF(line_char(i)/="")THEN
		        READ(line_char(i),*) bcVal(jj,j,nk)
		    ELSE
		        bcVal(jj,j,nk) = 0.
		    END IF

		    i = 4; j = 4
		    IF(line_char(i)/="")THEN
		        READ(line_char(i),*) bcVal(jj,j,nk)
		    ELSE
		        bcVal(jj,j,nk) = 0.
		    END IF

		    i = 5; j = 1
		    IF(line_char(i)/="")THEN
		        READ(line_char(i),*) bcVal(jj,j,nk)
		    ELSE
		        bcVal(jj,j,nk) = 0.
		    END IF

		    i = 6; j = 5
		    IF(line_char(i)/="")THEN
		        READ(line_char(i),*) bcVal(jj,j,nk)
		    ELSE
		        bcVal(jj,j,nk) = 0.
		    END IF

		    i = 7; j = 6
		    IF(line_char(i)/="")THEN
		        READ(line_char(i),*) bcVal(jj,j,nk)
		    ELSE
		        bcVal(jj,j,nk) = 0.
		    END IF

		    i = 8; j = 7
		    IF(line_char(i)/="")THEN
		        READ(line_char(i),*) bcVal(jj,j,nk)
		    ELSE
		        bcVal(jj,j,nk) = 0.
		    END IF

		    i = 10; j = 8
		    IF(line_char(i)/="")THEN
		        READ(line_char(i),*) bcVal(jj,j,nk)
		    ELSE
		        bcVal(jj,j,nk) = 0.
		    END IF

		    i = 14; j = 9  ! # gates
		    IF(line_char(i)/="")THEN
		        READ(line_char(i),*) bcVal(jj,j,nk)
		    ELSE
		        bcVal(jj,j,nk) = 0.
		    END IF

		    bcVal(jj,2,nk)=bcVal(jj,2,nk)*bcVal(jj,9,nk)   ! combining all gates into one
		    bcDes(jj,nk)=xcDes(ic+1,1)                       ! remember ic is the u/s cross section now
            continue		
 		ENDIF
!    ELSEIF(xcType(ic)==6) THEN    !  lateral weir
!        write(*,*) "unknown type of data"
!        write(*,*) "please contact Victor Huang to add internal boundary conditions for lateral weir."
!		cycle main_loop
	CASE DEFAULT
!        write(*,*) "unknown type of data"
!        write(*,*) "please contact Victor Huang to add additional internal boundary conditions."
!		stop
		cycle main_loop
	END SELECT
ENDDO main_loop

icm=ic
nriv=jj

!**********************************************************************
! read some options
!**********************************************************************
!WRITE(*,'(a,$)') 'Do you want to include floodplains in the numerical simulation(Y/N)? '
!READ(*,*) answer
answer = 'y'
IF(answer(1:1)=='Y' .or. answer(1:1)=='y') THEN
    floodPlain=.true.  
ELSEIF(answer(1:1)=='N' .or. answer(1:1)=='n') THEN
    floodPlain=.false.  
ELSE
    write(*,*) "wrong answer, program terminated!"
	stop
ENDIF

WRITE(*,'(a,$)') 'Do you want to skip some cross sections (Y/N)? '
READ(*,*) answer
!answer = 'n'
IF(answer(1:1)=='Y' .or. answer(1:1)=='y') THEN
    skipxs=.true.  
ELSEIF(answer(1:1)=='N' .or. answer(1:1)=='n') THEN
    skipxs=.false.  
ELSE
    write(*,*) "wrong answer, program terminated!"
	stop
ENDIF


IF(skipxs) THEN
	WRITE(*,'(a,$)') 'The geometry will be output in every N cross sections, input N: '
	READ(*,*) icInc 
	IF(icInc<=0) THEN
		write(*,*) "wrong answer, program terminated!"
		stop
	ENDIF   
ENDIF

IF(icInc>1) THEN
	ikept=0
	keep: DO
		IF(ikept==0) WRITE(*,'(a,$)') 'Do you want to keep a spectific cross sections (Y/N)? '
		IF(ikept>0) WRITE(*,'(a,$)') 'Do you want to keep another spectific cross sections (Y/N)? '
		READ(*,*) answer
		IF(answer(1:1)=='Y' .or. answer(1:1)=='y') THEN
			keepXS=.true.  
		ELSEIF(answer(1:1)=='N' .or. answer(1:1)=='n') THEN
			keepXS=.false.  
			EXIT keep
		ELSE
			write(*,*) "wrong answer, program terminated!"
			stop
		ENDIF

		IF(keepXs) THEN
			ikept=ikept+1
			WRITE(*,'(a,$)') 'The geometry of cross section M will be output, input M: '
			READ(*,*) xsKept(ikept) 
			IF(xsKept(ikept)<1 .or. xsKept(ikept)>icm) THEN
				write(*,*) "out of range, program terminated!"
				stop
			ENDIF   
		ENDIF
	ENDDO keep
ENDIF


!WRITE(*,'(a,$)') 'Do you want to read gis file (Y/N)? '
!READ(*,*) answer
answer = 'n'
IF(answer(1:1)=='Y' .or. answer(1:1)=='y') THEN
    read_gisfile = .true.  
ELSEIF(answer(1:1)=='N' .or. answer(1:1)=='n') THEN
    read_gisfile = .false.  
ELSE
    write(*,*) "wrong answer, program terminated!"
	stop
ENDIF

!read_gisfile = .true.
IF(read_gisfile)THEN
	DO I=1,10
		WRITE (*,'("Enter the file containing GIS info: ",$)')
		READ (*,'(A)') gisfile
		gisfile = ADJUSTL(gisfile)
		INQUIRE (file=gisfile,EXIST=exists)
		IF (.NOT.exists) THEN
			WRITE (*,'(2a/)') '  <<ERROR>> FILE NOT FOUND: ', gisfile
			CYCLE
			IF(i==10)THEN
				WRITE (*,'(2a/)') ' EXITING....'
				STOP
			END IF
		ELSE
			EXIT
		END IF
	ENDDO
	lugis = 21
	OPEN(UNIT=lugis,file=gisfile,status='old',IOSTAT=istat,ACTION='READ')
	IF(istat/=0)THEN
		WRITE(*,*) 'ERROR opening file',infile 
		STOP
	END IF

	! read header line
	CALL commchk(line,lineno,ndata,lugis,luo1)
	! read data
	read_gis:DO 
		CALL commchk(line,lineno,ndata,lugis,luo1)
		IF(line(1:6) == 'END OF')THEN
			EXIT read_gis
		END IF
		READ(LINE,*) itemp1,xcname,x1,y1,z1,des
		read(xcname,*,IOSTAT=istat)rm2
! find xc in HEC-RAS data
		RIVER0: DO jj=1,nriv
		    STATION0: DO ic=istart(jj),iend(jj)
				sdum = xcid(ic)
				read(sdum,*,IOSTAT=istat)rm1
				IF(sdum(1:7)==xcname(1:7) .or. ABS(rm1-rm2)<0.0002)THEN
					j = n_gis(ic)+1
					n_gis(ic) = j
					x_gis(j,ic) = x1
					y_gis(j,ic) = y1
				END IF
			END DO STATION0
		END DO RIVER0
	END DO read_gis
	CLOSE(lugis)
END IF

WRITE(*,'(a,$)') 'Do you wannt to set a maximum distance between points in a cross section (Y/N)? '
READ(*,*) answer
IF(answer(1:1)=='Y' .or. answer(1:1)=='y') THEN
	addpoints=.true.
	WRITE(*,'(a,$)') 'Please input the maximum horizontal distance between points: '
	READ(*,*) ymax 
	IF(ymax<0.5) THEN
		write(*,*) "Maxinum distance too small, program terminated!"
		stop
	ENDIF   
ELSEIF(answer(1:1)=='N' .or. answer(1:1)=='n') THEN
    addpoints=.false.  
ELSE
    write(*,*) "wrong answer, program terminated!"
	stop
ENDIF
! add points
! addpoints = .false.
IF(addpoints)THEN
    DO ic=1,icm 
        ! icm is the total number of cross sections 
        call insertPints(npt(ic),maxpts,ymax,y(:,ic),z(:,ic))
    ENDDO
END IF

WRITE(*,'(a,$)') 'Do you want to lower points if you have a flat bottom at some cross sections (Y/N)? '
READ(*,*) answer
!answer = 'n'
IF(answer(1:1)=='Y' .or. answer(1:1)=='y') THEN
    addLowerPoint=.true.  
ELSEIF(answer(1:1)=='N' .or. answer(1:1)=='n') THEN
    addLowerPoint=.false.  
ELSE
    write(*,*) "wrong answer, program terminated!"
	stop
ENDIF

deltZ=0.5
ymax=50.0
IF(addLowerPoint)THEN
    DO ic=1,icm 
        ! icm is the total number of cross sections 
        call insertLowerPoints(npt(ic),maxpts,bankL(ic),bankR(ic),ymax,deltZ,y(:,ic),z(:,ic))
    ENDDO
END IF



!**********************************************************************
! if floodplains are not needed, eliminate them
!**********************************************************************
DO IC=1,ICM
	bank1: DO j=1,npt(ic)
		IF(y(j,ic)>=bankl(ic)-TINY) THEN
			jbankL(ic)=j
			EXIT bank1
		ENDIF 
	ENDDO bank1
	bank2: DO j=npt(ic),1,-1
		IF(y(j,ic)<=bankr(ic)+TINY) THEN
			jBankR(ic)=j
			EXIT bank2
		ENDIF 
	ENDDO bank2
ENDDO


IF(floodPlain) THEN
    js=1
	je=npt
ELSE
    js=jbankL
	je=jbankR
ENDIF


!**********************************************************************
! output SRH-1D file
!**********************************************************************
RIVER1: DO jj=1,nriv
!   if there is a cross section associated with a internal boundary, keep it
    STATION1: DO ic=istart(jj),iend(jj)
!		DO nk=1,nki(jj)
!			IF(ic==nxi(jj,nk)) THEN
!		        icOutput(jj)=icOutput(jj)+1
!		        outXS(ic)=icOutput(jj)
!				CYCLE STATION1
!			ENDIF
!		ENDDO 
		IF(npt(ic)<=0) THEN
            outXS(ic)=0
		ELSEIF(MOD(ic-istart(jj),icInc)==0 .or. ic==icm) THEN
			icOutput(jj)=icOutput(jj)+1
			outXS(ic)=icOutput(jj)
		ELSEIF(ikept>0)THEN
			DO i=1,ikept
				IF(ic==xsKept(i)) THEN
			        icOutput(jj)=icOutput(jj)+1
			        outXS(ic)=icOutput(jj)
				ENDIF
			ENDDO 
		ELSE
            outXS(ic)=0
		ENDIF
	ENDDO STATION1
ENDDO RIVER1

min_dis = 0.
WRITE(luo2,'(100a)') ('*', j=1,100)
WRITE(luo2,'(100a)') 'Begin Internal Sturctures'
WRITE(luo2,'(100a)') ('*', j=1,100)

RIVER3: DO jj=1,nriv
    WRITE(luo2,'(a)')'********************************************'
    WRITE(luo2,'(a,2(a10,1x),i10)')'***	',rivDes(jj),reaDes(jj),jj
    WRITE(luo2,'(a3,1x,i10,1x,a)')'NAM',jj,reaDes(jj)
	WRITE(luo2,'(3h***,1x,5(A14,1x))') 'NKI(J)','b.c.','for','internal','station'
	WRITE(luo2,'(3hINF,1x,i14,1x)') nki(jj)
	WRITE(luo2,*)
    DO nk=1,nki(jj)
	    IF(ki(jj,nk)==8) THEN    ! radial gate
			WRITE(luo2,'(3h***,1x,A)') bcDes(jj,nk)
			WRITE(luo2,'(3h***,1x,4(A14,1x),A10,a17,I2)') 'NXI(J,nk)','KI(J,nk)','XTI(J,NK)','**********',bcName(jj,nk),'of internal B.C. ',nk
			WRITE(luo2,'(3hIFB,1x,2(I14,1x),F14.2)') outXS(nxi(jj,nk)),ki(jj,nk),xti(jj,nk)
			WRITE(luo2,'(3h***,1x,10(A14,1x))') 'C','W','T','Zsp','TE','BE','HE','Cw','GDIR','GTYPE'
			WRITE(luo2,'(3hI08,1x,8(F14.2,1x),2(I14,1x))') bcVal(jj,1:8,nk),0,0		    
			WRITE(luo2,'(3h***,1x,7(A14,1x))') 'WSEOpen','WSEClose','OpenRate','CloseRate','MaxOpen','MinOpen','InitOpen'
			WRITE(luo2,'(3hI8B,1x,A,1x)') '    please input gate opening function here'
	        WRITE(luo2,*)			
		ENDIF
	ENDDO
END DO RIVER3 

WRITE(luo2,'(100a)') ('*', j=1,100)
WRITE(luo2,'(100a)') 'Begin Cross Section Geometry'
WRITE(luo2,'(100a)') ('*', j=1,100)
RIVER4: DO jj=1,nriv
	WRITE(luo2,'(100a)') ('*', j=1,100)
    WRITE(luo2,'(a,2(a10,1x),i10)')'***	Begin River: ',rivDes(jj),reaDes(jj),jj
    WRITE(luo2,'(a3,1x,i10,1x,a)')'NAM',jj,reaDes(jj)

	!	skip cross section that channel length is 0
    DO ic=istart(jj),iend(jj)-1
		IF(x(ic)<min_dis)THEN
			outXS(ic+1) = 0
			xl(ic) = xl(ic) + xl(ic+1)
			xr(ic) = xr(ic) + xr(ic+1)
			x(ic) =  x(ic) +  x(ic+1)
		END IF
	END DO

 	!	add channel lenth to upstream cross section if this cross section is skipped
   DO ic=iend(jj),istart(jj),-1
		IF(outXS(ic) == 0)THEN
			xl(ic-1) =  xl(ic-1) + xl(ic)
			xr(ic-1) =  xr(ic-1) + xr(ic)
			x (ic-1) =  x (ic-1) + x (ic)
		END IF
	END DO



!   output cross sections
    STATION4: DO ic=istart(jj),iend(jj)
        IF(outXS(ic)==0) THEN
			cycle STATION4
		END IF
		WRITE(luo2,'(3h***,a)') '**********************************************************************************************************************************************************'
		WRITE(luo2,'(3h***,1x,a,a14,a)') 'station: ',xcid(ic)
		WRITE(luo2,'(3h***,1x,a)') xcDes(ic,1)
		WRITE(luo2,'(3h***,1x,a)') xcDes(ic,2)
		WRITE(luo2,'(3h***,1x,3(A,1x,I4,1x))') 'cross section',outXS(ic), 'of', icOutput(jj), 'at original',ic
		WRITE(luo2,'(3h***,1x,3(A14,1x))') 'XCID','ZDI','QDI'
		WRITE(luo2,'(3hXIN,1x,a30,1x,10(F14.2,1x))') xcid(ic),0.0, 0.0   !wse(ic),123.0
		WRITE(luo2,'(3h***,1x,7(A14,1x))') 'location','ch_l','ch_r','bec','ninterp','iHotC','base_rate'
		xloc = x(iend(jj))+xt(iend(jj))-xt(ic)
		!IF(jj>1)xloc = xloc + xt(iend(jj-1))

		WRITE(luo2,'(3hXST,1x,3(f14.2,1x),10(i14,1x))') xloc,xl(ic),xr(ic),0,0,0.,0
		
		WRITE(luo2,'(3h***,1x,3(A14,1x))') 'station','elevation','data'
		WRITE(luo2,'((3hXSP,1x,10(f14.2,1x)))') (y(j,ic),z(j,ic),j=js(ic),je(ic))

		IF(ndead(ic)>0) THEN
			WRITE(luo2,'(3h***,1x,3(A14,1x))') 'locl_dead','locr_dead','hdead'
		    WRITE(luo2,'((3hXIX,1x,9(f14.2,1x)))') (deadl(j,ic),deadr(j,ic),hdead(j,ic),j=1,ndead(ic))
		ENDIF
		IF(npdead(ic)>0) THEN
			WRITE(luo2,'(3h***,1x,3(A14,1x))') 'locl_pdead','locr_pdead','hpdead'
			WRITE(luo2,'((3hXPX,1x,9(f14.2,1x)))') (pdeadl(j,ic),pdeadr(j,ic),hpdead(j,ic),j=1,npdead(ic))
		ENDIF
		IF(nlev(ic)>0) THEN
			WRITE(luo2,'(3h***,1x,3(A14,1x))') 'locl_lev','locr_lev','hlev'
			WRITE(luo2,'((3hXLX,1x,9(f14.2,1x)))') (leveel(j,ic),leveer(j,ic),hlev(j,ic),j=1,nlev(ic))
		ENDIF
		IF(nblock(ic)>0) THEN
			WRITE(luo2,'(3h***,1x,3(A14,1x))') 'locl_block','locr_block','hblock'
		    WRITE(luo2,'((3hXBX,1x,9(f14.2,1x)))') (blockl(j,ic),blockr(j,ic),hblock(j,ic),j=1,nblock(ic))
		ENDIF

		WRITE(luo2,'(3h***,1x,2(A14,1x))') 'xloc_rcoef',	'rcoef'
		WRITE(luo2,'((3hXRH,1x,10(f14.3,1x)))') (loc_rr(j,ic),sect_rr(j,ic),j=1,num_rr(ic))
		WRITE(luo2,'(3h***,1x,2(A14,1x))') 'bankl',	'bankr'
		WRITE(luo2,'(3hXOX,1x,2(f14.2,1x))') bankl(ic),bankr(ic)
		IF(skew(ic)>TINY)THEN
			WRITE(luo2,'(3h***,1x,2(A14,1x))') 'skew'
			WRITE(luo2,'(3hXSK,1x,2(f14.2,1x))') skew(ic)
		END IF
		WRITE(luo2,'(3h***,1x,3(A14,1x))') 'expansion','contraction','minor'
		WRITE(luo2,'(3hXFL,1x,3(f14.2,1x))') fkec1(ic),fkec2(ic),fminor(ic)
		WRITE(luo2,'(3h***,1x,4(A14,1x))') 'GIS Cut Line'
!		IF(n_gis(ic)<1)THEN
!			n_gis(ic) = n_gis(ic-1)
!			DO j =1,n_gis(ic)
!				x_gis(j,ic) = x_gis(j,ic-1)
!				y_gis(j,ic) = y_gis(j,ic-1)
!			END DO
!		ENDIF
		IF(n_gis(ic)>1)THEN
			WRITE(luo2,'((3hXSL,1x,4(f14.2,1x)))') (x_gis(j,ic),y_gis(j,ic),j=1,n_gis(ic))
		ELSE
			WRITE(luo2,'((3hXSL,1x,4(f14.2,1x)))') y(1,ic),xloc,y(npt(ic),ic),xloc
		END IF
    ENDDO STATION4
ENDDO RIVER4


!**********************************************************************
! output thalweg
!**********************************************************************
RIVER5: DO jj=1,nriv
!   output cross sections
    STATION5: DO ic=istart(jj),iend(jj)
		if(outXS(ic)==0) cycle STATION5
		WRITE(luo4,'(i4,6F14.2,2x,2a8,3a40)') ic, xt(iend(jj))-xt(ic),MINVAL(z(jbankL(ic):jbankR(ic),ic)), &
&                                             x_gis(1,ic),y_gis(1,ic),x_gis(n_gis(ic),ic),y_gis(n_gis(ic),ic),  &            
&                                             "station:",xcDes(ic,1), xcid(ic),rivDes(jj),reaDes(jj)
    ENDDO STATION5
ENDDO RIVER5
write(*,*) "Program exculted successfully"

END

!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
SUBROUTINE commchk(line,lineno,ndata,lu,luo)
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
IMPLICIT NONE
!-----------------------------------------------------------------------------!
!                                                                             !
!  This subroutine reads the input file line by line.  The behavior is the    !
!  following:                                                                 !
!  - if the first character is a blank, skip this line and read next line;    !
!  - if the first character is a '*', print the line to the .DBG file and     !
!    read the next line;                                                      !
!  - if the first two characters are 'CM' (or 'cm'), print the line to the    !
!    .DBG file and read the next line (same behavior as '*');                 !
!  - if the first two lines are 'LW' or 'LX' (or 'lw' or 'lx'), ignore the    !
!    line and read the next line.                                             !
!                                                                             !
!-----------------------------------------------------------------------------!
! Dummy arguments.
INTEGER, INTENT(INOUT) :: lineno ! Line number of line read, for debugging.
CHARACTER(LEN=300), INTENT(OUT) :: line ! Line read from file.
CHARACTER TAB
DATA TAB/9/
INTEGER, INTENT(IN) :: lu,luo   ! input and output files, respectively
INTEGER, INTENT(out) :: ndata  ! number of data


CHARACTER, EXTERNAL :: to_upper
integer :: i,istat

line = ''
lineno = lineno + 1
READ (lu,'(A)',IOSTAT=istat) line
IF(istat/=0)THEN  ! error in read or end of file
!	WRITE(*,*) 'EOF encountered'
	line = 'END OF FILE'
	RETURN
END IF

DO i=1,300
	line(i:i) = to_upper(line(i:i))   ! Convert record name to upper case.
ENDDO
! Comment lines are printed out to .out file:
ndata=0
DO i=2,300
     if(line(i:i) == '\0') exit
     if((line(i-1:i-1) ==' ' .or. line(i-1:i-1) ==TAB .or. line(i-1:i-1)==',')    &
&       .and. (line(i:i) /=' ' .and. line(i:i) /=TAB .and. line(i:i)/=','    )) ndata=ndata+1
ENDDO

RETURN

END SUBROUTINE commchk

!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
CHARACTER FUNCTION to_upper(ch)
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  IMPLICIT NONE

! Converts a single character from lower case to upper case.  Leaves character
! unchanged if the character is not a lower case letter.
! F. Simoes, Dec. 2001

  ! Dummy argument.
  CHARACTER, INTENT(IN) :: ch

  ! Local variables.
  INTEGER, PARAMETER :: strand = IACHAR("A") - IACHAR("a")

  IF ("a" <= ch .AND. ch <= "z") THEN
    ! Lower case: change to upper case.
    to_upper = ACHAR(IACHAR(ch) + strand)
  ELSE
    ! Do nothing.
    to_upper = ch
  END IF

END FUNCTION to_upper

!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
SUBROUTINE data_err(id,line,lineno,luo)
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
IMPLICIT NONE
  INTEGER, INTENT(IN) :: luo   ! output files, respectively

! Print-out of error message when a wrong record is found and abort program.

! Dummy arguments.
  INTEGER, INTENT(IN) :: lineno
  CHARACTER(LEN=*), INTENT(IN) :: id
  CHARACTER(LEN=300), INTENT(IN) :: line

  WRITE (luo,'(2X,"Expecting record ",A," - Found:",/,2X,A,/,2X,"in line", &
        I6,/,"Please recheck the dataset.")') id,line,lineno
  WRITE (*,'(2X,"Expecting record ",A," - Found:",/,2X,A,/,2X,"in line", &
        I6,/,"Please recheck the dataset.")') id,line,lineno
	PAUSE
  STOP

END SUBROUTINE data_err


!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
SUBROUTINE spaceLine(num,luo)
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
IMPLICIT NONE
INTEGER, INTENT(IN) :: num,luo   ! output files, respectively

INTEGER :: i
DO i=1,num
	WRITE(luo,*)
ENDDO
END SUBROUTINE spaceLine




